module cev_finder_optimizer_mod
USE functions_mod
USE parameters_mod
IMPLICIT NONE

CONTAINS



!!!!!TOTAL

SUBROUTINE cev_finder(cev_out,sim_consumption_start,sim_energy_start,sim_labor_start,survive,&
    sim_consumption_end,sim_energy_end,sim_labor_end,numberN,non_energy)
INTEGER, intent(in)   :: numberN,non_energy
DOUBLE PRECISION,    intent(in)   :: sim_consumption_start(J,numberN),sim_energy_start(J,numberN)
DOUBLE PRECISION,    intent(in)   :: sim_labor_start(J,numberN),survive(J)
DOUBLE PRECISION,    intent(in)   :: sim_consumption_end(J,numberN),sim_energy_end(J,numberN)
DOUBLE PRECISION,    intent(in)   :: sim_labor_end(J,numberN)
DOUBLE PRECISION,    intent(out)  :: cev_out
DOUBLE PRECISION :: diff_out



diff_out=brent(-50.0D0, .5D0, 50.0D0,cev_differ,.000000001D0,cev_out,&
    sim_consumption_start,sim_energy_start,sim_labor_start,survive,&
    sim_consumption_end,sim_energy_end,sim_labor_end,numberN,non_energy)

cev_out=-cev_out
end subroutine cev_finder


!This is objective function to find lambda0
DOUBLE PRECISION FUNCTION cev_differ(cev_test,sim_consumption_start,sim_energy_start,sim_labor_start,survive,&
    sim_consumption_end,sim_energy_end,sim_labor_end,numberN,non_energy)
INTEGER, intent(in)   :: numberN,non_energy
    DOUBLE PRECISION,    intent(in)   :: sim_consumption_start(J,numberN),sim_energy_start(J,numberN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor_start(J,numberN),survive(J)
    DOUBLE PRECISION,    intent(in)   :: sim_consumption_end(J,numberN),sim_energy_end(J,numberN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor_end(J,numberN),cev_test



    DOUBLE PRECISION	::value_start(J,numberN),value_end(J,numberN)

    INTEGER ::it,it2
    !! Taxes

    do it =J,1,-1
        do it2=1,numberN
            if (it==J) then
                value_end(it,it2)=utility_last_cev(sim_consumption_end(it,it2),sim_energy_end(it,it2),adult_equivalence(it),cev_test,non_energy)
                value_start(it,it2)=utility_last_cev(sim_consumption_start(it,it2),sim_energy_start(it,it2),adult_equivalence(it),0.0D0,non_energy)

            elseif(it>Jnr-1) then
                value_end(it,it2)=valuer_ret_cev(sim_consumption_end(it,it2),sim_energy_end(it,it2),survive(it),value_end(it+1,it2),adult_equivalence(it),cev_test,non_energy)
                value_start(it,it2)=valuer_ret_cev(sim_consumption_start(it,it2),sim_energy_start(it,it2),survive(it),value_start(it+1,it2),adult_equivalence(it),0.0D0,non_energy)


            else
                value_end(it,it2)=valuer_cev(sim_consumption_end(it,it2),sim_energy_end(it,it2),sim_labor_end(it,it2),survive(it),value_end(it+1,it2),adult_equivalence(it),cev_test,non_energy)
                value_start(it,it2)=valuer_cev(sim_consumption_start(it,it2),sim_energy_start(it,it2),sim_labor_start(it,it2),survive(it),value_start(it+1,it2),adult_equivalence(it),0.0D0,non_energy)


            end if
        end do
    end do
    cev_differ=abs(sum(value_end(1,:))-sum(value_start(1,:)))
end FUNCTION cev_differ




DOUBLE PRECISION FUNCTION brent(ax,bx,cx,func,tol,xmin,sim_consumption_start,sim_energy_start,sim_labor_start,survive,&
    sim_consumption_end,sim_energy_end,sim_labor_end,numberN,non_energy)
USE nrtype; USE nrutil, ONLY : nrerror
IMPLICIT NONE
INTEGER, intent(in)   :: numberN,non_energy
    DOUBLE PRECISION, INTENT(IN)    :: ax,bx,cx,tol
    DOUBLE PRECISION,    intent(in)   :: sim_consumption_start(J,numberN),sim_energy_start(J,numberN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor_start(J,numberN),survive(J)
    DOUBLE PRECISION,    intent(in)   :: sim_consumption_end(J,numberN),sim_energy_end(J,numberN)
    DOUBLE PRECISION,    intent(in)   :: sim_labor_end(J,numberN)
    DOUBLE PRECISION, INTENT(OUT)   :: xmin
INTEGER, PARAMETER          :: ITMAX=500
DOUBLE PRECISION, PARAMETER ::CGOLD=0.3819660D0,ZEPS=1.0e-3_sp*epsilon(ax)
INTEGER                     ::iter
DOUBLE PRECISION            :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
DOUBLE PRECISION :: func
a=min(ax,cx)
b=max(ax,cx)
v=bx
w=v
x=v
e=0.0
fx=func(x,sim_consumption_start,sim_energy_start,sim_labor_start,survive,&
    sim_consumption_end,sim_energy_end,sim_labor_end,numberN,non_energy)
fv=fx
fw=fx
do iter=1,ITMAX
    xm=0.5D0*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0D0*tol1
    if (abs(x-xm) <= (tol2-0.5D0*(b-a))) then
        xmin=x
        brent=fx
        RETURN
    end if
    if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.0D0*(q-r)
        if (q > 0.0D0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5D0*q*etemp) .or. &
            p <= q*(a-x) .or. p >= q*(b-x)) then
            e=merge(a-x,b-x, x >= xm )
            d=CGOLD*e
        else
            d=p/q
            u=x+d
            if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
    else
        e=merge(a-x,b-x, x >= xm )
        d=CGOLD*e
    end if
    u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
    fu=func(u,sim_consumption_start,sim_energy_start,sim_labor_start,survive,&
        sim_consumption_end,sim_energy_end,sim_labor_end,numberN,non_energy)
    if (fu <= fx) then
        if (u >= x) then
            a=x
        else
            b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
    else
        if (u < x) then
            a=u
        else
            b=u
        end if
        if (fu <= fw .or. w == x) then
            v=w
            fv=fw
            w=u
            fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
            v=u
            fv=fu
        end if
    end if
end do
call nrerror('brent lambda: exceed maximum iterations')
CONTAINS
!BL
SUBROUTINE shft(a,b,c,d)
    DOUBLE PRECISION, INTENT(OUT)   :: a
    DOUBLE PRECISION, INTENT(INOUT) :: b,c
    DOUBLE PRECISION, INTENT(IN)    :: d
a=b
b=c
c=d
END SUBROUTINE shft
END FUNCTION brent












end module cev_finder_optimizer_mod
